In [1]:
from functools import partial
from pathlib import Path
In [2]:
import pandas as pd
import numpy as np
In [3]:
import plotly.offline as py
py.init_notebook_mode()
In [4]:
import holoviews as hv
hv.extension('plotly')
In [5]:
import matplotlib.pyplot as plt
In [6]:
%matplotlib inline
In [7]:
# first value is default
# rest are alternative value for each realization
PARAMS = {
    'STELLAR_BARYON_FRAC': (0.05, 0.001, 1.),
    'STELLAR_BARYON_PL': (0.5, -0.5, 1.),
    'ESC_FRAC': (0.01, 0.001, 1.),
    'ESC_PL': (-0.5, -1., 0.5),
    'M_TURNOVER': (5e8, 1e8, 1e10),
    't_STAR': (0.5, 0., 1.),
    'L_X': (40.5, 38., 42.)
}
In [8]:
from dautil.plot import iplot_column_slider, plot_column_slider
In [9]:
def read_global_evolution(path):
    df = pd.read_csv(path, delim_whitespace=True, header=None, names=('zp', 'filling_factor_of_HI_zp', 'Tk_ave', 'x_e_ave', 'Ts_ave', 'T_cmb*(1+zp)', 'J_alpha_ave', 'xalpha_ave', 'Xheat_ave', 'Xion_ave'))
    df.set_index('zp', inplace=True)
    return df
In [10]:
def parse_case(string):
    for level, idx in enumerate(map(int, string.split('_'))):
        if idx != 0:
            break
    if idx == 0:
        return 'default'
    else:
        param = list(PARAMS)[level]
        return '{}={:.6}'.format(param, PARAMS[param][idx])
In [11]:
basedir = Path('~/21cmfast/21cmFAST/Parameter_spaces').expanduser()
In [28]:
df_path = pd.DataFrame((basedir.glob('**/global*')), columns=['path'])
In [30]:
df_path['case'] = df_path.path.map(lambda path: path.parent.parent.parent.name)
In [32]:
df_path['name'] = df_path.case.map(parse_case)
In [35]:
# filter out bad cases
df_path = df_path[~df_path.case.isin(('0_0_0_0_0_1_0',))]
In [36]:
df_path
Out[36]:
path case name
0 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_0_0_0_0 default
1 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 1_0_0_0_0_0_0 STELLAR_BARYON_FRAC=0.001
2 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 2_0_0_0_0_0_0 STELLAR_BARYON_FRAC=1.0
3 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_1_0_0_0_0_0 STELLAR_BARYON_PL=-0.5
4 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_2_0_0_0_0_0 STELLAR_BARYON_PL=1.0
5 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_1_0_0_0_0 ESC_FRAC=0.001
6 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_2_0_0_0_0 ESC_FRAC=1.0
7 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_1_0_0_0 ESC_PL=-1.0
8 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_2_0_0_0 ESC_PL=0.5
9 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_0_1_0_0 M_TURNOVER=1e+08
10 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_0_2_0_0 M_TURNOVER=1e+10
12 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_0_0_2_0 t_STAR=1.0
13 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_0_0_0_1 L_X=38.0
14 /home/kolen/21cmfast/21cmFAST/Parameter_spaces... 0_0_0_0_0_0_2 L_X=42.0
In [65]:
df = pd.concat(
    (read_global_evolution(path).x_e_ave for path in df_path.path),
    axis=1,
#     keys=df_path.name
    keys=df_path.case
)
In [39]:
df.sort_index(inplace=True)
In [ ]:
df.to_hdf(
    '21cmfast-x_e_ave.hdf5',
    'df',
    format='table',
    complevel=9,
    fletcher32=True
)
df = pd.read_hdf('21cmfast-x_e_ave.hdf5')
In [13]:
py.iplot(iplot_column_slider(df))
In [14]:
py.plot(iplot_column_slider(df), filename='x_e.html')
Out[14]:
'file:///Users/kolen/git/source/UCB-PHYS229-2019Spring-final-project/notebook/x_e.html'
In [53]:
df.head()
Out[53]:
name default STELLAR_BARYON_FRAC=0.001 STELLAR_BARYON_FRAC=1.0 STELLAR_BARYON_PL=-0.5 STELLAR_BARYON_PL=1.0 ESC_FRAC=0.001 ESC_FRAC=1.0 ESC_PL=-1.0 ESC_PL=0.5 M_TURNOVER=1e+08 M_TURNOVER=1e+10 t_STAR=1.0 L_X=38.0 L_X=42.0
zp
6.000593 0.000767 0.000222 0.009026 0.002381 0.001066 0.000737 0.001194 0.000811 0.000811 0.000893 0.000477 0.000492 0.000175 0.011135
6.140605 0.000706 0.000219 0.008592 0.002244 0.000951 0.000681 0.001104 0.000744 0.000737 0.000823 0.000441 0.000461 0.000174 0.010294
6.283417 0.000651 0.000215 0.008146 0.002103 0.000850 0.000630 0.001020 0.000683 0.000672 0.000759 0.000410 0.000433 0.000172 0.009487
6.429086 0.000601 0.000212 0.007688 0.001959 0.000761 0.000583 0.000941 0.000628 0.000615 0.000701 0.000382 0.000408 0.000171 0.008719
6.577667 0.000556 0.000209 0.007220 0.001816 0.000682 0.000541 0.000868 0.000579 0.000565 0.000648 0.000358 0.000385 0.000170 0.007983

would be z values when x_e crosses 1 if interpolated linearly

In [62]:
x_0 = df.index[0]
dx = df.index[1] - x_0
y_0 = df.iloc[0]
dy = df.iloc[1] - y_0
(1. - y_0) * dx / dy + x_0
Out[62]:
name
default                      -2287.171042
STELLAR_BARYON_FRAC=0.001   -46040.369232
STELLAR_BARYON_FRAC=1.0       -313.397730
STELLAR_BARYON_PL=-0.5       -1013.149663
STELLAR_BARYON_PL=1.0        -1212.503555
ESC_FRAC=0.001               -2496.161195
ESC_FRAC=1.0                 -1552.367862
ESC_PL=-1.0                  -2079.405430
ESC_PL=0.5                   -1873.745363
M_TURNOVER=1e+08             -1988.517394
M_TURNOVER=1e+10             -3922.750146
t_STAR=1.0                   -4490.596661
L_X=38.0                    -86331.398001
L_X=42.0                      -158.671414
dtype: float64

Generate class .ini files

In [117]:
INI = '''H0 = 67.32117
omega_b = 0.02238280
N_ur = 2.03066666667
omega_cdm = 0.1201075
N_ncdm = 1
omega_ncdm = 0.0006451439
YHe = 0.2454006
tau_reio = 0.05430842
n_s = 0.9660499
A_s = 2.100549e-09
non linear = halofit
output = tCl,pCl,lCl,mPk
lensing = yes
root = output/{case}-
write warnings = yes
write parameters = yes
input_verbose = 1
background_verbose = 1
thermodynamics_verbose = 1
perturbations_verbose = 1
transfer_verbose = 1
primordial_verbose = 1
spectra_verbose = 1
nonlinear_verbose = 1
lensing_verbose = 1
output_verbose = 1
format = camb
reio_parametrization = reio_inter
reio_inter_num = {reio_inter_num}
reio_inter_z = {reio_inter_z}
reio_inter_xe = {reio_inter_xe}
write thermodynamics = yes
'''

Binned

In [66]:
def binning(array, binwidth=4):
    n, m = array.shape
    n_trunc = (n // binwidth) * binwidth
    if n_trunc != n:
        print(f'array truncated from {n} to {n_trunc}.')
        array = array[:n_trunc]
    return array.T.reshape(m, -1, binwidth).mean(axis=-1)
In [105]:
def format_float_array_trunc(array):
    string = ','.join(f'{i:.15}' for i in array)
    print(len(string))
    return string
In [83]:
def format_float_array_st(array):
    string = ','.join(f'{i:.16E}' for i in array)
    print(len(string))
    return string
In [81]:
def format_float_array(array):
    string = ','.join(map(str, array))
    print(len(string))
    return string
zp, x_e = binning(df_temp[['zp', 'default']].values, binwidth=10)zp.shape[0]print(format_float_array(zp))print(format_float_array(x_e))hv.Curve(np.column_stack((zp, x_e)))

Unbinned

In [67]:
# setting boundary values according to documentation on reio_inter
df.loc[0.] = -2.
df.loc[6.] = -1.
df.loc[30.] = 0.
In [69]:
df.sort_index(inplace=True)
In [71]:
df.head()
Out[71]:
case 0_0_0_0_0_0_0 1_0_0_0_0_0_0 2_0_0_0_0_0_0 0_1_0_0_0_0_0 0_2_0_0_0_0_0 0_0_1_0_0_0_0 0_0_2_0_0_0_0 0_0_0_1_0_0_0 0_0_0_2_0_0_0 0_0_0_0_1_0_0 0_0_0_0_2_0_0 0_0_0_0_0_2_0 0_0_0_0_0_0_1 0_0_0_0_0_0_2
zp
0.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000 -2.000000
6.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000 -1.000000
6.000593 0.000767 0.000222 0.009026 0.002381 0.001066 0.000737 0.001194 0.000811 0.000811 0.000893 0.000477 0.000492 0.000175 0.011135
6.140605 0.000706 0.000219 0.008592 0.002244 0.000951 0.000681 0.001104 0.000744 0.000737 0.000823 0.000441 0.000461 0.000174 0.010294
6.283417 0.000651 0.000215 0.008146 0.002103 0.000850 0.000630 0.001020 0.000683 0.000672 0.000759 0.000410 0.000433 0.000172 0.009487
In [82]:
 
850
Out[82]:
'0.0,6.0,6.000593,6.140605,6.283417,6.429086,6.577667,6.72922,6.883805000000001,7.041481,7.202311,7.366357000000001,7.533683,7.704358,7.878445,8.056014,8.237134,8.421877,8.610314,8.802521,8.998571,9.198542999999999,9.402514,9.610563,9.822775,10.03923,10.260015,10.485215,10.714919,10.949218,11.188202,11.431966000000001,11.680605,11.934217,12.192902,12.456759,12.725895,13.000413,13.280420999999999,13.56603,13.85735,14.154497000000001,14.457587,14.766739000000001,15.082072,15.403714,15.731788,16.066423,16.407751,16.755905,17.111023,17.473244,17.842709,18.219563,18.603954,18.996032999999997,19.395954,19.803873,20.21995,20.644348,21.077234,21.51878,21.969154,22.428537,22.897108,23.37505,23.862551,24.359802,24.866999,25.384338,25.912025,26.450266,26.999271000000004,27.559258000000003,28.130444,28.713053000000002,29.307314,29.913459999999997,30.0'
In [74]:
for name, series in df.items():
    break
In [84]:
series.head()
Out[84]:
zp
0.000000   -2.000000
6.000000   -1.000000
6.000593    0.000767
6.140605    0.000706
6.283417    0.000651
Name: 0_0_0_0_0_0_0, dtype: float64
In [106]:
for name, series in df.items():
    assert len(format_float_array_trunc(series.values)) < 1024
994
995
973
980
988
993
990
992
993
991
990
993
991
967
In [123]:
outbasedir= Path('~/21cmfast/21cmFAST/Parameter_spaces').expanduser()
In [124]:
outbasedir
Out[124]:
PosixPath('/home/kolen/21cmfast/21cmFAST/Parameter_spaces')
In [125]:
reio_inter_num = df.shape[0]
reio_inter_z = format_float_array(df.index.values)
for case, series in df.items():
    reio_inter_xe = format_float_array_trunc(series.values)
    with open(outbasedir / f'{case}.ini', 'w') as f:
        print(INI.format(case=case, reio_inter_num=reio_inter_num, reio_inter_z=reio_inter_z, reio_inter_xe=reio_inter_xe), file=f)
850
994
995
973
980
988
993
990
992
993
991
990
993
991
967
In [122]:
for case in df:
    (outbasedir / f'{case}.ini').unlink()